unit ppot;

interface
uses common;

{Global declarations for pair potentials}
Const
   NDTA	= T_ELE; {chemical types; comes from common.p}
   NDTT	= NDTA*(NDTA+1) div 2; {num_of_pairs}
   NDVR	= 2000; {potentials grid dimension}
   nele	= z_ele; {110, number of elements}
Type
   fr	= array[1..ndvr] of myreal;	
Var
   RR		    : fr;
   vr,vr1,vr2,vr3   : array[1..ndtt] of fr;
   rcut,frcut,rcut2 : myreal;
   rcplus,cueff	    : myreal;
   hr,rstart,epot   : myreal;
   mtt		    : array[1..nele,1..nele] of integer;
   pab		    : array[1..ndtt,1..2] of integer;
   ipotp	    : array[0..nele] of integer;
   chmap1	    : integer;
   iext		    : array[1..3] of integer;
   nvr		    : integer;

procedure upp(i1,i2 :integer; ar : myreal; var v,vd,vdd:myreal);
procedure uppv(i1,i2 :integer; ar : myreal; var v:myreal);
procedure DSPLIN(X,Y : fr; N:integer; B1,BN:myreal; Var Y1,Y2,Y3:fr);
procedure POTIN;
procedure pp_inp;
procedure unit_cell(multi :integer);

implementation

procedure unit_cell(multi :integer);
var txt	: text;
   i,j	: integer;
begin
   if debug[10] then write('Reading INPUT/PPOT/uc.inp: ');
   reset(txt,'INPUT/PPOT/uc.inp');
   readln(txt,lpar[1],lpar[2]);
   close(txt);
   if debug[10] then writeln(lpar[1]:1:4,' ',lpar[2]:1:4);
   for i:=1 to 3 do
      for j:=1 to 3 do bl3[i,j]:=0.0;
   for i:=1 to 2 do
      for j:=1 to 2 do bl3[i,j]:=ba[i,j]*lpar[1];
   lpar[2]:=lpar[2]*multi;
   bl3[3,3]:=lpar[2];
   for i:=1 to 3 do {transposed}
      for j:=1 to 3 do blt[i,j]:=bl3[j,i];
end;
procedure upp(i1,i2 :integer; ar : myreal; var v,vd,vdd:myreal);
{note: hr is global}
var ir,ijt : integer;
   dr	   : myreal;
begin
   ijt:=mtt[i1,i2];
   if (ijt>0) and (ijt<=ndtt) then begin
      if (ar<rstart) then begin
	 v:=1e3; vd:=0.0; vdd:=0.0;
	 {writeln('upp-short: ',i1:1,'-',i2:1,'   ',ar:1:4);}
      end else begin
	 IR:=trunc((AR-rstart)/HR)+1;
	 DR:=(AR-RR[IR]);
	 V:=VR[ijt,ir]+VR1[ijt,ir]*DR+VR2[ijt,ir]*DR*DR+VR3[ijt,ir]*DR*DR*DR;
	 VD:=((3.0*VR3[ijt,ir]*DR+2.0*VR2[ijt,ir])*DR+VR1[ijt,ir])/AR;
	 VDD:=6.0*VR3[ijt,ir]*DR+2.0*VR2[ijt,ir];
	 {if ar<2.3 then
	 writeln(ir:5,ar:8:4,dr:8:4,i1:5,i2:3,ijt:5,vr[ijt,ir]:12:7,v:12:7);}
      end;
   end else begin
      write('UPP: mismatch atom-type PP-type, ');
      writeln('ijt out of range: ',i1:3,i2:3,ijt:1);
      v:=0.0; vd:=0.0; vdd:=0.0;
      halt;
   end;
end; { upp }   
procedure uppv(i1,i2 : integer; ar : myreal; var v:myreal);
{returns only energy, no derivatives}
{note: hr is global}
var
   ir,ijt : integer;
   dr	   : myreal;
begin
   ijt:=mtt[i1,i2];
   if (ijt>0) and (ijt<=ndtt) then begin
      if (ar<rstart) then begin
	 v:=1e3;
	 {writeln('upp-short: ',i1:1,'-',i2:1,'   ',ar:1:4);}
      end else begin
	 IR:=trunc((AR-rstart)/HR)+1;
	 DR:=(AR-RR[IR]);
	 V:=VR[ijt,ir]+VR1[ijt,ir]*DR+VR2[ijt,ir]*DR*DR+VR3[ijt,ir]*DR*DR*DR;
	 {if ar<2.3 then
	 writeln(ir:5,ar:8:4,dr:8:4,i1:5,i2:3,ijt:5,vr[ijt,ir]:12:7,v:12:7);}
      end;
   end else begin
      write('UPP: mismatch atom-type PP-type, ');
      writeln('ijt out of range: ',i1:3,i2:3,ijt:1);
      v:=0.0;
      halt;
   end;
end; { upp }   
   
procedure DSPLIN(X,Y : fr; N:integer; B1,BN:myreal; Var Y1,Y2,Y3:fr); 
{
[AFTER MK, fortran]:
 C     PROGRAM FOR CALCULATING COEFICIENTS FOR EVALUATING FUNCTION GIVEN
 C     FOR A SET OF POINTS BY CUBIC SPLINE 
 C     X - ARRAY OF POINTS    
 C     Y - ARRAY OF FUNCTION VALUES 
 C     N - NUMBER OF POINTS 
 C     B1,BN - SECOND DERIVATIVE OF THE FUNCTION AT THE END POINTS
 C     SET B1:=BN:=0 IF NO INFORMATION IS  AVAILABLE
 C     Y1,Y2,Y3 - ARRAYS OF SPLINE COEFICIENTS 
 C     VALUE OF FUNCTION AT ARBITRARY  X  CAN BE EVALUTED BY FORMULA
 C     F[X]:=[[[Y3[K]*[X-X[K]]+Y2[K]]*[X-X[K]]+Y1[K]]*[X-X[K]+Y[K]
 C     WHERE K IS THAT OF MAXIMAL X[K] LESS THAN X
 }
Var
   m2,i,j : integer;
   r,s    : myreal;
begin
   M2:=N-1;
   S:=0.0;
   for I:=1 to M2 do begin
      Y1[I]:=X[I+1]-X[I];
      R:=(Y[I+1]-Y[I])/Y1[I];
      Y2[I]:=R-S;
      S:=R;
   end;
   S:=0.0;
   R:=0.0;
   Y2[1]:=B1;
   Y2[N]:=BN; 
   for I:=2 to M2 do begin
      Y2[I]:=Y2[I]+R*Y2[I-1];
      Y3[I]:=2.0*(X[I-1]-X[I+1])-R*S;
      S:=Y1[I];
      R:=S/Y3[I];
   end;
   for J:=2 to M2 do begin
      I:=M2+2-J;
      Y2[I]:=(Y1[I]*Y2[I+1]-Y2[I])/Y3[I];
   end;
   for I:=1 to M2 do begin
      S:=Y1[I];
      R:=Y2[I+1]-Y2[I];
      Y3[I]:=R/S;
      Y2[I]:=3.0*Y2[I];
      Y1[I]:=(Y[I+1]-Y[I])/S-(Y2[I]+R)*S;
   end;
end;
procedure POTIN;
{READS AND TRANSFORMS POTENTIAL TABLES}
Var
   ir,io,jo,o,i,j    : integer;
   dr,scr,sce,dr0,r0 : myreal;
   yc,fi,y,z,b1,bn   : myreal;
   ar,v,vd,vdd,eps   : myreal;
   rnew		     : myreal;
   it,npp,noneq	     : integer;
   txt		     : text;
   wr1,wr2,wr3	     : fr;
   vtmp		     : array[1..ndtt] of fr;

   procedure smooth_and_spline;
   Var ir,it : integer;
   begin           
      R0:=FRCUT*RCUT;
      YC:=RCUT/R0;
      for IR:=1 to NVR do begin {smooth cutoff}
	 FI:=1.0;
	 Y:=RR[IR]/R0;
	 Z:=(Y-YC)/(1.0-YC);
	 FI:=3.0*Z*Z*Z*Z-8.0*Z*Z*Z+6.0*Z*Z;
	 IF (Y<1.0+eps) THEN FI:=1.0;
	 IF (Y>YC+eps) THEN FI:=0.0;
	 for IT:=1 to npp do VR[it,ir]:=VR[it,ir]*FI;
      end;
      for IT:=1 to npp do begin
	 VR[it,1]:=VR[it,2]*2.0-VR[it,3];
	 b1:=0.0; bn:=0.0;
	 DSPLIN(RR,VR[it],NVR,b1,bn,vr1[it],vr2[it],vr3[it]);
      end;
   end; { smooth_and_spline }			     
begin
   noneq:=0;
   eps:=1e-4;
   for i:=1 to nele do
      for j:=1 to nele do begin
	 mtt[i,j]:=0;
      end;
   if debug[10] then write('Reading INPUT/PPOT/pp.tab: ');
   reset(txt,'INPUT/PPOT/pp.tab');
   readln(txt,nvr,scr,npp,sce);
   if ndtt<npp then begin
      writeln('Insufficient PP dimension NDTA.');
      halt;
   end;
   for o:=1 to npp do begin
      read(txt,i,j); mtt[i,j]:=o; mtt[j,i]:=o;
      pab[o,1]:=i; pab[o,2]:=j;
   end; readln(txt);
   for IR:=1 to NVR do begin
      READ(txt,rr[ir]);
      rr[ir]:=rr[ir]*scr;
      if (ir>1) then dr0:=rr[ir]-rr[ir-1];
      if (ir=2) then dr:=dr0;
      if (ir>2) then if (abs(dr-dr0)>1e-5) then begin
	 {writeln('dr(1) dr(',ir:1,') : ',dr:1:5,' ',dr0:1:5);
	 writeln('PP table: non-equidistant mesh not implemented.');
	 halt;}
	 noneq:=1; {nonequidistant mesh}
      end;
      for it:=1 to npp do begin
	 READ(txt,VR[it,ir]); 
	 vr[it,ir]:=vr[it,ir]*sce;
      end; readln(txt);
   end;
   close(txt);
   rstart:=rr[1];
   SMOOTH_AND_SPLINE;
   HR:=DR; {will be modifd for noneq mesh}
   {make equidistant mesh if the current is nonequidistant}
   if(noneq=1)then begin
      {we take say 2-ce as much points as defined}
      ir:=1; {pointer to noneq r value}
      hr:=(rr[nvr]-rr[1])/(5*nvr);
      for i:=0 to 5*nvr do begin
	 rnew:=rstart+i*hr;
	 dr:=rnew-rr[ir];
	 if((ir+1<nvr)and(dr>rr[ir+1]-rr[ir]))then begin
	    ir:=ir+1;
	    dr:=rnew-rr[ir];
	 end;
	 for o:=1 to npp do begin
	    vtmp[o,i+1]:=VR[o,ir]+VR1[o,ir]*DR+VR2[o,ir]*DR*DR+VR3[o,ir]*DR*DR*DR;
	 end;
	 {writeln(i:1,' ',ir:1,' ',dr:1:4,' ',rnew:1:5,' ',vtmp[1,i+1]:1:6);}
      end;
      nvr:=5*nvr+1;
      for i:=1 to nvr do
	 for o:=1 to npp do vr[o,i]:=vtmp[o,i];
      for i:=1 to nvr do rr[i]:=rstart+(i-1)*hr;
      SMOOTH_AND_SPLINE;
   end;
   {rcut:=rcut*YC;}
   rewrite(txt,'INPUT/pp.out');
   dr0:=hr/5.0;
   for ir:=0 to 5*nvr do begin
      ar:=rstart+dr0*ir;
      write(txt,ar:1:3,' ');
      for o:=1 to npp do begin
	 upp(pab[o,1],pab[o,2],ar,v,vd,vdd);
	 write(txt,v:18,' ');
      end; writeln(txt);
   end;
   close(txt);
   if (debug[10]) then writeln(nvr:1,' r-points.');
end;
   
procedure pp_inp;
var txt	: text;
var i,o	: integer;
begin
   for i:=0 to nele do ipotp[i]:=0;
   if (debug[10]) then write('Reading INPUT/PPOT/pp.inp: ');
   reset(txt,'INPUT/PPOT/pp.inp');
   for o:=1 to 3 do read(txt,iext[o]); readln(txt);
   readln(txt,rcut); rcut2:=rcut*rcut;
   readln(txt,frcut);
   close(txt);
   if debug[10] then writeln('rcut=',rcut:1:2,', frcut=',frcut:1:3);
   potin;
end; { pp_inp }

end.
